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Abstract. We review the black hole solutions of the ghost-free massive gravity 
theory and its bimetric extension and outline the main results on the stability of 
these solutions against small perturbations. Massive (bi)-gravity accommodates exact 
black hole solutions, analogous to those of General Relativity. In addition to these 
solutions, hairy black holes - solutions with no correspondent in General Relativity - 
have been found numerically, whose existence is a natural consequence of the absence 
of Birkhoff’s theorem in these theories. The existence of extra propagating degrees of 
freedom, makes the stability properties of these black holes richer and more complex 
than those of General Relativity. In particular, the bi-Schwarzschild black hole exhibits 
an unstable spherically symmetric mode, while the bi-Kerr geometry is also generically 
unstable, both against the spherical mode and against superradiant instabilities. If 
astrophysical black holes are described by these solutions, the superradiant instability 
of the Kerr solution imposes stringent bounds on the graviton mass. 
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Massive gravity is a modification of General Relativity (GR) based on the idea of 
equipping the graviton with mass. A model of non self-interacting massive gravitons 
was first suggested by Fierz and Pauli in the early beginnings of field theory 0 0. 
This original model was subsequently shown by van Dam, Veltman and Zakharov 
(vDVZ) to differ from GR even at small distance scales, ruling out the theory on the 
basis of Solar system tests [3]. A solution to this problem was later conjectured by 
Vainshtein [3], who argued that GR could be recovered at small distances by including 
non-linear terms in the field equations of the hypothetical massive gravity theory. 
Much later, rigorous studies of several non-linear completions of massive gravity showed 
that this would indeed generically be the case (for a review see Ref. [5j). However, 
generic nonlinear versions of the Fierz-Pauli theory, although able to recover GR via 
the Vainshtein mechanism, turned out to reveal another pathology — the so-called 
Boulware-Deser ghost |6j. This ghost problem has only been solved recently in a series 
of works 0 0 El EDI El El fl3j . where it was shown that for a subclass of massive 
potentials the Boulware-Deser ghost does not appear, both in massive gravity — a 
theory with one dynamical and one fixed metric, the so-called de Rham, Gabadadze 
and Tolley (dRGT) model — and its bi-gravity extension (see Refs. [13, [E] for recent 
reviews on massive gravity||J The bi-gravity theory contains two dynamical metrics, 
which interact with each other via non-derivative terms. If ordinary matter only couples 
to one of the metrics, then one can interpret the theory as an extension of Einstein’s 
gravity with an extra spin-2 field coupled to gravity in a particular non-minimal way. 

Despite its very recent formulation, a considerable amount of work has been done 
to understand the implications of the ghost-free theory. Although a great deal of the 
interest has been focused on building cosmological models out of these theories (see [15] 
and references therein), an increasing effort is underway to construct and understand 
the properties of black hole (BH) solutions in massive (bi)-gravity. Understanding BHs 
in any theory of gravity is obviously of extreme importance, not only to understand the 
highly non-linear regime of the theory, but also to possibly look for deviations from GR 
(see e.g. Ref. [IS])- 

It is interesting that the first BH solutions in a nonlinear theory of massive gravity 
were found in the context of high-energy physics mi. with no mention on the possible 
application of the theory as a modification of gravity. Reference mi introduced non¬ 
bidiagonal solutions of a bi-metric theory, with both metrics having the Schwarzschild- 
(anti)de-Sitter form, but not diagonal simultaneously. However, BHs in bi-gravity did 
not attract much attention until the discovery of the massive (bi)-gravity theory safe 

f A settled term often used in the literature to refer to dRGT theory and its extensions is the “ghost-free 
massive gravity”. This name should be used with care, since dRGT theory is free from the Boulware- 
Deser ghost (or Ostrogradski ghost), but the other degrees of freedom are not necessary healthy for 
particular solutions. Therefore the right term should rather be “Boulware-Deser ghost free massive 
gravity” or “Ostrogradski ghost free massive gravity”, although we will also sometimes refer to dRGT 
theory and its extension simply as “ghost free massive gravity”, keeping in mind the above reservation. 


Black holes in massive gravity 


3 


from the Boulware-Deser ghost. The discovery of this theory led to a plethora of works on 
the subject, and similar solutions have been put forward for both dRGT massive gravity 
and its bi-gravity extension mmmmm- More interestingly, it has been shown that 
spherically symmetric solutions which do not exist in GR are also present [20-, [23]. These 
are the only examples known so far of “hairy” solutions in theories of massive gravity 
(see Refs. [23, 25] for reviews covering this topic). Later, Ref. |26] generalized some of 
these previous works finding a class of non-bidiagonal Reissner-Nordstrdm solutions in 
dRGT massive (bi)-gravity, while Ref. [27] was able to construct a family of rotating 
BHs in the same theories. 

Along with the search for BH solutions in massive (bi)-gravity, in the last couple 
of years progress has been made to understand the stability properties of these 
solutions. One of the most striking results of this ongoing effort was the conclusion 
that the bidiagonal Schwarzschild solution is unstable, as found in [28] and confirmed 
independently in |29j . In contrast, the same type of instability has been shown to be 
absent for non-bidiagonal solutions Pi. 

As a by-product of studying the stability of these solutions, by understanding 
how small fluctuations behave in theories with massive gravitons it is also possible to 
understand how gravitational waveforms might differ from GR. The extra gravitational 
polarizations and the nontrivial dispersion introduced by a putative small graviton 
mass may leave important imprints on gravitational waveforms from, e.g., inspiralling 
compact objects. Understanding these effects is thus necessary given that advanced 
gravitational-wave detectors mm will begin operation (Advanced LIGO [31] is in 
fact expected to begin to collect data in mid-2015) and the first direct detection of 
gravitational waves is expected to take place within the next decade. With this is mind, 
Refs. [29, .33] studied generic perturbations of Schwarzschild and Kerr BHs in massive 
(bi)-gra.vity, with the particularly interesting result that the Kerr solution would be 
prone to another kind of instability, related to the superradiant scattering of bosonic 
fields with spinning BHs (see Refs. [33, 35] for recent reviews). This instability, which 
was known to occur for massive scalar [351 EZ] and vector [38, GES 3D] perturbations of 
rotating BHs in GR, was shown to occur in massive gravity due to the natural existence 
of linear massive spin-2 perturbations in the theory. If BHs in the Universe were to 
be described by these solutions, the most striking consequence of this instability is the 
existence of a bound in the graviton mass [29]. 

This work is meant to review in some detail BH solutions in dRGT massive gravity 
and its bi-gravity extension, including both exact and numerical solutions, with the 
focus on the stability issues — a fresh and fast developing topic. We introduce in Sec. [2] 
massive gravity theories: Fierz-Pauli theory, non-linear completions and the dRGT 
model with its bi-gravity extension. Then, in Sec. [3] we present recent developments 
on BH solutions in dRGT model: the exact bidiagonal and non-bidiagonal solutions, 
as well as numerical solutions. Sec. [4] is devoted to results of the past two years on 
perturbations and stability of BHs in massive (bi)-gravity. In Sec. [5] we conclude with 
open issues and work in progress. 
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2.1. Fierz- Pauli theory 


Before introducing the theory of massive (bi)-gravity, we briefly review the linear massive 
gravity model, also known as the Fierz-Panli theory 0 El- The linear theory can 
also be viewed as an expansion of the full non-linear massive gravity model around a 
Minkowski background. Expanding the Einstein-Hilbert action in metric perturbations 
as g fW = r + h pV) where is the Minkowski metric rj )MJ = diag(-l, +1, +1, +1) 
and the indices of h pv are moved up and down with the metric g pu , and keeping only 
quadratic terms in the action we obtain the linear approximation of GR, 


Sgr — 




( 1 ) 


where Mp is the Planck mass and 

+ ]^d p d p h p + ^ dpdyh’l - ^ry lv (d p d' J h p(T - Dh), 

is the linearized Einstein tensor G pu = £ lw + 0(h 3 ). When matter is present, the 
metric perturbation h pv is also coupled to the energy-momentum tensor T pu , via the 
interaction term h^T^, but since here we will mostly consider vacuum solutions, we 
omit this term. The action (|TJ) contains only derivative terms. By adding non-derivative 
quadratic terms h 2 (where h = h ILU rf lIJ ) and h pv h pv to the action (JT|) , a linear massive 
gravity theory is obtained. If the non-derivative terms are added in a special combination 
oc {h liv h ,iu — h 2 ), the action takes the Fierz-Panli form p, 2], 
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where m is the mass parameter, corresponding to the graviton mass. This particular 
mass combination describes the only consistent linear Lorentz-invariant theory for a 
massive spin-2 field. Models with other mass terms can be shown to necessarily contain 
a physical ghost degree of freedom (DEI. which is in turn related to the Ostrogradski 
ghost HD. The Fierz-Panli theory ([2]) fails to pass the Solar system tests of gravity, due 
to the fact that a massive graviton has more polarizations than a massles one, therefore 
modifying the gravitational interaction. Even in the zero graviton mass limit this leads 
to a difference between the massless and massive gravitation interaction which is known 
as the vDVZ discontinuity [3]. 


2.2. Non-linear massive gravity 

Starting from the quadratic action ([2]), one can try to guess a non-linear generalization 
of the theory, in the same way that full GR is a non-linear generalization of the quadratic 
action (JT|) . Obviously, the first term in (J2]) constructed out of the metric g should be 
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the Einstein-Hilbert term[§} However, it is not possible to get a non-linear massive term, 
using only the metric g^, since the only nontrivial term corresponds to a Lagrangian 
density proportional ^—g, which stands for the cosmological constant. Thus, one way to 
have a non-trivial mass term is to add a second metric, say / M „, which can be chosen to 
be fixed (in this case the theory has a preferred background, i.e. “aether”) or dynamical, 
in which case the theory is called bi-gravity or massive bimetric gravity. The metric g liv 
can be non-derivatively coupled to the second metric f /w . in order to form a non-trivial 
mass term. Non-linear mass terms should be chosen such that: the action is invariant 
under a coordinate change common to both metrics; there is a (almost) flat solution for 
g tJ „] in the limit where g tiV = g^u + h^ and = rg w the potential at quadratic order for 
h^ v takes the Fierz-Pauli form (|2]). In spite of the restrictions formulated above, there is 
a huge freedom in choosing the mass term. In fact one can choose the interaction term 
in a class of functions satisfying these conditions (see e.g. Ref. H). The term, 

V=9 (a.» - - Ur) (<r <T - s'Vh. 

considered in [45], is an example of such an interaction. As it has been proposed 
by Vainshtein [4], the inclusion of non-linear terms in the equations of motion for 
a massive graviton help to solve the problem of the vDVZ discontinuity. However, 
as was later noticed in Ref. [6] a proper proof of this conjecture was lacking. The 
status of the Vainshtein mechanism in non-linear massive gravity was in fact only 
solved recently, when analytic and numerical studies confirmed that the Vainshtein 
mechanism indeed works in this model [46] (for a recent review see [5]). Furthermore 
in the work by Boulware and Deser [6] another serious problem has been found: the 
non-linear interaction terms generically lead to ghost instabilities (appearing at the 
non-linear level). The Boulware-Deser instability can be thought as another face of the 
Ostrogradski ghost instability [ 41] . 

2.3. dRGT gravity 

Although, in general, the Boulware-Deser ghost persists in non-linear massive and bi¬ 
gravity theories, it has been found by de Rham, Gabadadze and Tolley (dRGT) that in 
the so-called “decoupling limit” — a limit where the degrees of freedom of the theory 
(almost) decouple — there is a restricted subclass of potential terms, for which the 
Boulware-Deser ghost is absent even at the non-linear level [3 U El E]. Later a 
full Hamiltonian analysis confirmed the absence of the Boulware-Deser ghost in this 
model HD1 EJ Q3], while fully covariant proofs were given in Ref. 03 for a subset of 
possible massive terms and for generic mass terms in [48J (see also the review [15]). 

To formulate the dRGT theory and its extension to the bi-gravity case, it is 
convenient to introduce functions e*, of matrices X , which represent the elementary 

§ In principle, since Lovelock’s theorem [42] is not valid in massive gravity, other kinetic terms could 
be possible. However it was shown in Refs. 1431 that any new kinetic term would generically lead to 
ghost instabilities. 
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symmetric polynomials of the eigenvalues of X. In the case of 4 x 4 matrices they are 
given by (cf. e.g. | 12 j), 


e„ = 1, e, = [A'], e 2 — t ([A ] 2 - [A 2 ]) , e 3 = i ([A ] 3 - 3[A][A 2 ] + 2[A 3 ]) , 
U = 4 ([A ] 4 - 6 [A] 2 [A 2 ] + 3[A 2 ] 2 + 8 [A][A 3 ] - 6 [A 4 ]) , 


(3) 


where [X] stands for the trace of A". Note that e 4 can also be written in the simpler 
way, e 4 = det(X). The building block of the dRGT mass term is the square root of the 
matrix g _1 f, i.e. one defines the matrix 7 as follows, 


7 = Vs *f, 


(4) 


with the matrix elements of 7 defined as 7 ^ = yWW- The action of the most general 
bi-gravity theory without the Boulware-Deser ghost then reads [49j . 


S = Mp 


J d 4 x^g 


k =4 


R[g\ - 2m 2 ^fte fc ( 7 ) 


k =0 


+ kM| / d 4 xV^Jn[f], 


(5) 


where R[g] and 7 Z[f] are the Ricci scalars for the metrics g and / respectively, /3 n and 
k are arbitrary coefficients. The bi-gravity action (J5]) contains both the kinetic term for 
the metric g and one for the metric /. The coefficient k marks the difference between 
the Planck masses for the g and / metrics. The original dRGT theory [2, El E] (dRGT 
massive gravity) is given by the action (|5]) dropping the last term, such that there is no 
dynamics for the metric /, contrary to the bi-gravity theory. Note also that the terms /3 0 
and /?4 do not give a mass to the graviton. Since a /—g e 4 ( 7 ) = y/—g det ( 7 ) = \/—f, the 
/3o term describes a cosmological constant for the metric g, while the /3 4 term corresponds 
to the cosmological constant of the metric /. Thus, there is a three parameter family of 
massive bi-gravity theories parametrized by /3k, with k — 1,2,3 (which becomes a two 
parameter family once the mass of the graviton is fixed). 

The action ([5]) can be written in an alternative form, using the matrix /C, defined 
as, 


K, = I- 7 . 


( 6 ) 


The bidiagonal extension of the dRGT action then reads, 

S bi G =M 2 P j d 4 x^g ( R[g} + 2m 2 U[g, /] + 2m 2 A g ) 

r (7) 

+ nM'p / d 4 x\f^-j (7 Z[f] + 2m 2 A/) , 

where 


U[g, f] = e 2 (K) + a 3 e 3 (K ) + a 4 e 4 (/l), 

with the following identifications for the coefficients, 

Po = — (A g + 6 + 4q?3 + a 4 ) , Pi — (3 + 307 + cr 4 ) , 

P2 — ~ (1 + 2a3 + a 4 ), P 3 = ( ot 3 + a 4 ), / 3 4 = — (kA f + a 4 ) , 


( 8 ) 
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and where in the action ([7]) we explicitly separated the cosmological terms. 

The equations of motion derived from (Jr]) by variation with respect to g^ v and / M „ 
read, 


G\ = m 2 (n - A s «) , 


G“„ = rn‘ 


-t K 




(9) 

( 10 ) 


where and Q ,L U are the Einstein tensors for the metrics g and / correspondingly and 
the energy-momentum tensors coming from the interaction terms have the form 

rp _ 7 , 0 SU 

-L liu m z 


5gi JV 

- <W ll (K - KK) + 039 ^ (w< - [K]K + (K 2 K) 
+ 049 ^la - IMC- + |/C](V 2 )“ - ( K.X) + U'Xu,- 


( 11 ) 


(K - [K«) - a 3 /,„ 7 : (u,j; - \k]k: + ac 2 )*) ( 12 ) 

- (%V - U'lK + [K](K 2 r - (AC 3 )S) . 


It can be shown that T^ u and are symmetric [58]. Note also a useful relation between 
the energy-momentum tensors, 


n = -T'‘„+M«. (13) 

where T ,l v is found from by raising an index with g^, while to get T^, one raises an 
index with the metric f^ v . Furthermore the Bianchi identity implies the conservation 
conditions 


V^ = 0, 


vjv = 0 , 


(14) 


where V ff and V f are the covariant derivatives with respect to g /iu and respectively. 
Note that due to Eq. (13) these two conditions are not independent. Finally, if one 
consider the metric / MI/ to be non-dynamical, then Eq. (10) must be excluded. 


3. Black holes in massive (bi)-gravity 

The structure of solutions in massive (bi)-gravity is more complex that in GR, mainly 
due to the fact that this theory has two metrics (see e.g. Ref. [50] to see how the global 
structure of these solutions is affected by the co-existence of two metrics). In particular, 
the well-known Birkhoff’s theorem for spherically symmetric solutions does not apply. 
This suggests that in massive (bi)-gravity the classes of BH solutions are richer than in 
GR, 

Indeed, the spherically symmetric BH solutions in bi-gravity theories can be divided 
into two classes. The first class corresponds to the case for which the metrics cannot be 
brought simultaneously to a bidiagonal form. Said differently, in this class, if one metric 
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in some coordinates is diagonal, the other metric is not. BHs of the second class have 
two metrics which can be both written in the diagonal form, but not necessarily equal. 

The first BH solutions for a nonlinear massive gravity theory (suffering from the 
Boulware-Deser ghost) were constructed in Ref. [17] . Much later, spherically symmetric 
solutions for other classes of ghosty massive bi-gravity theories were found and classified 
in detail in Refs. [501 EQ. In the framework of the original dRGT model a class of non¬ 
bidiagonal Schwarzschild-de-Sitter solutions was found in Ref. [TS] . In Refs. |[T5l [20]. 
spherically symmetric BH solutions were found for the bi-metric extension of dRGT 
theory, while spherically symmetric (charged and uncharged) solutions in the dRGT 
model for a special choice of the parameters of the action were presented in Refs. J2T, '25|. 
More recently, Ref. [26] found a more general class of charged BH solutions (in both the 
dRGT model and its bi-gravity extension). Finally Ref. [27J generalized these findings 
by including rotation in the geometry. This last class of solutions, jointly with the 
charged solutions of Ref. [26J, includes as particular cases most of the previously found 
spherically symmetric solution^]] 

Interestingly, spherically symmetric BH solutions with hair — solutions differing 
from the Schwarzschild family — were also found in bi-gravity theory, both with Anti-de 
Sitter and flat asymptotics [23]. 

Some good reviews on solutions of BHs in massive gravity already exist, e.g., 
Refs. [23, [25]. Here we will mainly focus in giving a different and unified treatment 
of all the solutions found so far. 

3.1. Analytic solutions in massive (bi)-gravity 

Due to the complexity of the field equations, it turns out that it is easier to find analytical 
BH solutions of the first class (with non-bidiagonal metrics). In this section we mostly 
focus on this type of BHs, although some bidiagonal solutions will also be presented. 
Since the BH solutions of the original dRGT model can easily be obtained from the 
bi-metric ones, all our calculations will be done for the bi-metric theory described by 
the action (J7]). 

To warm up, let us demonstrate the simplest BH solutions. We notice that for 
9fiu = fpuv the potential term in the action ([7]) is zero, U[g, f] =0. Comparison with GR 
then guarantee that in vacuum this theory admits the Kerr-(Anti) de Sitter metric as 
a solution. Taking for simplicity A g = Af = 0, and considering the particular case of 

| With the exception of Schwarzschild non-bidiagonal solutions, presented in [T5] . where an extra 
constant of integration appears in the solution. However, in [ 501 it has been argued (for similar BH 
solutions in a ghostly massive gravity) that the extra parameter should be set to a specific value for 
the solutions to be physical. In this case the solutions of m are a particular subclass of the solutions 
found in In 125], a method was presented to construct more general spherically symmetric non¬ 

bidiagonal solutions. These solutions are implicitly written in terms of one function (of two coordinates), 
which must satisfy a particular PDE, and thus describe an infinite-dimensional family of solutions (a 
similar technique has been used in [52] to find de Sitter solutions in dRGT massive gravity). 
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static solutions one then finds that the two metrics with the same Schwarzschild form, 

dr 2 


ds 2 g = ds) = 


is a 


(i - /) ^ + 

solution of the bi-gravity theory (JtJ) . 


1 - r " 


+ r W, 


(15) 


3.1.1. Static spherically symmetric solutions. To find other solutions it proves to be 
much more convenient to use coordinates which are regular at the horizon. One reason 
is that using such coordinates, we automatically avoid problems related to the singular 
behavior of the Schwarzschild coordinates at the horizon. In fact, for bidiagonal solutions 
this becomes a physical singularity, except for the particular case where the two metrics 


share the same horizon, like in the solution given in Eq. (15). More on this subject can 
be found in Ref. [53] (see also Ref. [53] where it was shown that the temperature of each 
horizon must also be the same). 

We start therefore with the advanced bi-Eddington-Finkelstein coordinates and 
assume the following ansatz for the two metrics 

1 — — — dv 2 + 2 dvdr + r 2 dQ 2 , (16) 


ds 2 g = 


ds) = C 2 


| 1 — — — I d yZ + 2 dvdr + r z dkl z 


(17) 


where r g and rf are the Schwarzschild radii for the g -and /-metrics correspondingly, l g 
and If are the “cosmological” radii and C is a constant. The Einstein tensors have the 
same form as in GR, 


G* = - S M (T 1 = 

^ v 12 VI =7 V 

L g 


CH) 


5v 


(18) 


while the energy-momentum tensor (fTT| takes the simple form, 


/ A, 


m m — 

V 


\ 


rpr 

1 v 

0 

0 


0 

A, 

0 

0 


0 

0 


0 \ 
0 
0 

As ) 


(19) 


where 


X g = —{C - 1 ) (0{C - l ) 2 - 3 a(C - 1 ) + 3) , 
is the effective cosmological constant for the metric g, 

T\ — —y (0(C — l ) 2 — 2a(C — 1) + 1) 

is the only non-diagonal term and we introduced the following constants, 
Q. = 1 “t~ 0 ( 3 , (3 = 0(3 0(4 . 



( 20 ) 


( 21 ) 


( 22 ) 


The energy-momentum tensor for the metric / can then easily be found from Eqs. (19) 


and (13). Clearly it must have the same form as Eq. (19), but with different values 
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for its components. Indeed, the diagonal part of consists of the components with 
values, 

Xf — C(C — 1) ((C - 1) 2 (1 - a + P) + 3(C - 1)(1 - a) + 3) , (23) 

while the off-diagonal term reads, 

T r v = -T\ . (24) 

To satisfy these modified Einstein’s equations with the effective energy-momentum 
tensor given in Eq. (19), the non-bidiagonal terms T r v and T T V must vanish. Depending 
on the way we “kill” those terms, different branches of solutions exist. We either need 
to put to zero the expression in the first or in the second parentheses of Eq. (21). First 
of all, one can easily see that the solution (15) is recovered for r g = rf, l g = If —$■ oo 
and C — 1 (in different coordinates though), since in this case = T IJ u = 0. A slightly 
more general solution is the bi-Schwarzschild-de-Sitter solution, which is obtained for 
r g = ry, l g = If, C = 1 and A g = Af — 3 /(ml g ) 2 . 

Another class of bidiagonal solutions arises for r g = r/, l g — If, but 1. In this 
case, the non-diagonal components of the energy-momentum tensors are zero, but the 
diagonal components X g and Xf are not, so from Eqs. i§. ( jToj ) we have the following 
relations, 

X f 


3 = (■ mlg ) (A g -Xg), 3 = C ( ml g y 


A/ C 4 K 


(25) 


For a given set of parameters in the action (|7j), these equations determine l g and C 
entering the metrics (16) and (17). 

The non-bidiagonal class of BH solutions emerges for rj ^ r g (and generically 
If 7 ^ lg). In order to “kill” the non-bidiagonal terms T r v and T r v we need to satisfy the 
relation, 

P(C - l) 2 - 2a(C - 1) + 1 = 0, (26) 

which fixes C in terms of the parameters of the action, while the de-Sitter radii are 
found from the field equations ([9]) and (10), which reduce to 

X f 


3 = (mlg) 2 (Ag - Xg), 3 = (^(m If) 2 ( A f 


C 4 k 


( 27 ) 


3.1.2. Black hole solutions with electric charge. Adding an electromagnetic source to 
the theory given by action (JtJ) , it is also possible to End charged BH solutions. We 
consider now the following action, 

S = S biG -^J d 4 x^Fg U F^, (28) 

where S b iG is given by (j7|) and the electromagnetic field is coupled to the metric g only. 
The field equation (|9]) for the ^-metric is then modified to include the energy-momentum 
tensor of the electromagnetic field, 

G\ = m 2 (n - A,K) + A (F m F„ a - 1stFX . (29) 
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while the field equation (10) for the metric / remains the same. We change accordingly 


the ansatz for the metric g as follows 


2 2 N 

ds 2 = — (l — — + -y- — yy) dv 2 + 2 dvdr + r 2 dQ J 


(30) 


while keeping the ansatz for the metric / ( [T7| ). 

With this modification the off-diagonal component of the energy-momentum 


tensor (21) is given by, 


n = -J b(C - l) 2 - 2 a(C - 1) + 1) ' r 


’ + '( 31 ) 


r 


Notice that the condition (26) still gives T r v = 0 as before, while the condition (27) 
provides the balance of the bare and the effective cosmological constants, fixing l g and 
If. Taking the vector potential in the standard form, 

A^iQ/r, 0,0,0}, (32) 

where Q is the charge of the g-BH, we find that the field equations are satisfied if the 
conditions ([26l). ([27b and 

(33) 


\/2Mpr q — Q , 


are also satisfied. To summarise, a class of charged BH solutions for the theory (28) 


is given by the metrics (30) and (17), with the vector potential (32) and the 


conditions (26), (27) and (33). 


3.1.3. Enhanced symmetry of black hole solutions. It turns out that for the particular 
combination of parameters of the action, f3 = a 2 , the space of solutions is wider than 
in the generic case fd ^ a 2 , due to an enhanced symmetry for this choice of parameter. 
Taking the following general ansatz for the metrics, 

ds 2 = g vv dv 2 + 2 g vr dvdr + g rr dr 2 + r 2 dQ 2 , (34) 

ds 2 = C 2 [f vv dv 2 + 2 f vr dvdr + f rr dr 2 + r 2 dQ 2 ] , (35) 

with g vv , g vr , g rr , f vv , f vr , f rr functions of v and r, with the conditions /3 = a 2 and 


(26), the energy-momentum tensors take the form of effective cosmological constants. 


Note that in the metrics (34) and (35) we only required spherical symmetry and not 


the precise ansatz for the metrics, in contrast to the Schwarzschild-de-Sitter ansatz (16) 


and ( [l7| ) in the study above. Of course, we still need to satisfy the modified Einstein’s 
equations with effective cosmological constants, thus the metrics g and / must have 
the Schwarzschild-de-Sitter form (or Reissner-Nordstrom-de Sitter, in case of non-zero 


charge) in some coordinates, as long as the conditions (27) are fulfilled. However, in 


this case, one of the metrics is not fixed with respect to the other, there is freedom to 
make a coordinate change of the form, 


v —> v(v, r), r —>• r(v, r), 


(36) 
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in one metric, without touching the other metric at the same time. One can perform 


an independent coordinate change (36) for each metric and the result is also a solution. 
Therefore, the solution for /3 = a 2 has an extra symmetry, which is absent in the general 
case. In the case with a fixed metric /, the fact that there is an extra freedom in this 
solution has been discussed in detail in Ref. j-TS . 

In fact, several solutions for BHs presented in the literature fall in this general class 
of solutions. In particular, the solutions of Ref. [22] and Ref. |21j . are examples of the 


solution described above, with particular coordinate changes of the form (36) for one of 
the metrics g and /, or both. More details can be found in Ref. 


3.1.4- Rotating black hole solutions The approach to look for BH solutions using the 
regular at the horizon metric ansatz also works — with some modifications — for 
rotating BHs [21] • Indeed, for the theory ([7]) let us assume the metrics g and / to 
have the form of the original Kerr metric element, 


ds 2 g = 


dsj = C 2 


where 


r q r\ , , 2 

-A- ( dv + a sin' 

P 2 J 

+2 {dv + a sin 2 6dq i) (dr + a 


sm' 


+ p 2 (dd 2 + 


- 1 - 


2 rfr 

P 2 

+2 (dv + a 


(dv + a sin 2 Odcf) ‘ 


sm 


(dr + a sin 2 ddq i) + p 2 (d6 2 + 


sm 


sm 


2 2,2 2 /i 

p = r + a cos 0, 


6d(j) 2 ) . 

(37) 

9d(f) 2 )~\ . 

(38) 


(39) 


and a is the rotation parameter. The ansatz (37) and (38) can be viewed as an extension 


of (16) and (|I7|) for rotating BHs. First of all, we note that since both metrics (37) and 


(38) coincide with the Kerr solution, then 

Gfiv — G/iv — 0 . 


(40) 


In our ansatz (37) and (38) the space for both metrics is asymptotically flat, in contrast 


to the previous cases (16), (17) and (30), where we allowed asymptotically (A)dS space- 


times. We have to set the right hand sides of (|9]) and (10) to zero in order to satisfy 
Einstein’s equations, although one could straightforwardly generalize these solutions to 
the asymptotically (A)dS case. 

Lengthy but straightforward calculation of the energy-momentum tensors gives 

0 0 0 \ 

{R)rp r 

0 


(R)rpfl — 


X, 

(R)rpr 


X, 


\ 


0 

0 


0 

0 


0 

Xn 


<t> 


0 


(41) 




where X g is given by (20) and 

{R) T\ = ~ (P(C - l) 2 - 2a(C - 1) + 1) (rg ~ r/)r \ 


p- 


( 42 ) 
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C 


(B) n = -17 (/S(C - l ) 2 - 2 a(C - 1) + 1) 


a(r g — Tf)r sin 2 9 


(43) 


2 v ' ' ■ ' p* 

Here we use the superscript (R) to emphasize that this corresponds to a rotating solution. 
Note that (42) coincides with (21) for a = 0 and l g — If — oo, while the r.li.s. of (43) is 
zero when there is no rotation, a = 0. From Eq. (13) it is now easy to find the 

off-diagonal components are given by, 


(R)'j-r ( R)rjnr _ _(. R)rjr, 




(44) 


while the diagonal components are equal to A f, given by Eq. (23). 

Both off-diagonal components of (41) are zero if the condition (26) is satisfied (as 
in the spherically symmetric case). On the other hand, the diagonal components, acting 
as an effective cosmological constant, must be balanced by an appropriately chosen A 9 
and A/, 

A / 


— A g, A f 


(45) 


which is simply the old condition (25) in the limit of asymptotically flat spacetimes, 
l g oo and If oo. 


We have just shown that the metrics (37) and (38) with C given by p 6 l) and A 




A f given by (45) are rotating BH solutions of the massive gravity model (|7j). 

It is not difficult now to get BH solutions of the dRGT model (with one fixed metric) 
from the BH solutions in bi-metric massive gravity. The dRGT action is given by g]) 
without the last term, which is the dynamical part of the second metric. Therefore, one 
simply needs to retract the equations of motion of the second metric, assuming each time 
that the second metric is flat [^} In particular, in the case of the Schwarzschild-de-Sitter 
BH, the solution of dRGT model is given by (16), and (17) with ry = 0 and If = oo, 
(so that fg, v is flat although written in the ingoing Eddington-Finkelstein coordinates), 
with the condition (26) and the Erst condition of (27). Similarly, a charged BH of the 
dRGT model is given by @, gg with Tf = 0 and If = oo, the condition ( [26] ) , and the 
first condition of (27) and, additionally, by (33). Finally, a rotating solution in dRGT is 
obtained from the solution of the bi-metric extension in a similar manner. In this case 


the metric / is given by (38) with r/ = 0 and If = oo. One can check that the /-metric is 
flat for this choice of parameters, although the line element is written in an unusual form. 
It can be obtained from the canonical Minkowski metric ds 2 M = — dt 2 + dx 2 + dy 2 + dz 2 
by the coordinate change, t — v — r, x + iy — (r — m)e*^sin 9, z = r cos#, and by the 
subsequent replacement r -)■ Cr, -)■ Cv, a —¥ Ca. In this case, the metric g is given 
by (37) where the condition (26) and the first constraint of (45) must be satisfied. 


One can also consider the non-dynamical metric to be non-flat, where each choice corresponds to a 
different theory. In the particular case where the non-dynamical metric is flat, bidiagonal BH solutions 
do not exist, since these are necessarily singular at the horizon [55] , 
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So far, all the solutions we showed belong to the Kerr family, which are also solutions of 
GR. However, the absence of no-hair theorems in massive gravity and in particular the 
non-validity of Birkhoff’s theorem for spherically symmetric spacetimes, suggests the 
existence of other solutions such as BHs endowed with massive spin-2 hair. Due to the 
complexity of the field equations (J9| and (10), finding such solutions requires the use of 
numerical methods. A detailed study of those solutions was performed in Ref. [20] who 
also studied the case where both metrics are diagonal but not necessarily proportional. 
In particular, asymptotically AdS hairy BH solutions were shown to exist. This was 
further extended in Ref. [23] where a family of asymptotically flat hairy BHs was found. 
More recently the same techniques were used to find vacuum wormhole solutions [56] in 
these theories. 

All the spherically symmetric uncharged solutions we discussed in the previous 
section are related through a coordinate change with the Schwarzschild-(A)dS metric. 
It is possible to find other solutions by considering static spherically symmetric solutions 
of the field equations ([9]) and (10), with the generic ansatz for the metrics given by 

g fll/ dx fl dx l/ = — Q 2 dt 2 + N~ 2 dr 2 + R 2 dfl 2 , (46) 


fuvdx^dx 1 ' = — a 2 dt 2 + b 2 dr 2 + U 2 dVt 2 


(47) 


where Q, N ,a ,b, R and U are radial functions. Gauge freedom allow us to 
reparametrize the radial coordinate r such that R(r ) = r. To simplify the field equations 
we also introduce the radial function Y(r) defined as b = U'/Y , where ' = d/dr. After 
using the conservation condition (14) the problem can be reduced to a closed system of 
three coupled first-order ODE’s for the functions N, Y and U (the derivation of these 
equations can be found in Ref. m and is also available online in a Mathematica 
notebook 0 ): 

N' = T\ (r, N, Y, U, /i, k, a 3 , aq) 

Y' = 2 - 2 (r, N, Y, U,/x,k, aq, aq) (48) 

U' = T 3 (r, N, Y, U, /a, k, aq, aq). 

The remaining two functions Q and a can then be evaluated from algebraic equations 
of the form: 


Q l Q' — R±(r, N,Y, U,h,k, a 3 , aq), (49) 

Q~ x a = Rr 0 (r,N,Y,U,ii,K,a 3 ,a A ). (50) 


Here we defined the parameter 



(51) 


which, as we will see later, corresponds to the graviton mass in some backgrounds. As 
discussed previously, when f^ v = C 2 g^ u , we find the same solutions of GR, namely 
Schwarzschild-(A)dS, where the value of the cosmological constant will depend on C 
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through Eq. (25). However, when we do not require the metrics to be proportional the 
complexity of the system of Eqs. (48) requires them to be solved numerically. 

To solve these equations, appropriate boundary conditions at the event horizon rh 
must be imposed. For the spacetime to be smooth at the horizon, bidiagonal solutions 
must share the same horizon [53], which implies that Q{rh) = N{rh ) = Y(rh) = ct{rh) = 
0. Assuming a power-series expansion at the horizon of the form 


n 2 = E' 


n> 1 


U = u r h + 


r h y 


c n \r 


V 2 = ^b n (r 


r h y 


n> 1 


n)’ 


(52) 

(53) 


n> 1 


it follows that, for a given set of parameters of the theory, the solutions are parametrized 
by one single free parameter u. Fixing the value of u = C, with C obtained from 
Eq. (25) and integrating numerically the equations from the horizon r = up to 
infinity, the solution is given by Schwarzschild-(A)dS, with the cosmological constant 
fixed by Eq. (25) [23] • On the other hand more interesting non-trivial solutions appear 


when choosing u = C + Su. 


3.2.1. Asymptotically AdS hairy solutions. Using this method hairy deformations of 
the Schwarzschild-AdS geometry were found in Ref. [20]. These solutions approach AdS 
spacetime when r —y oo, showing deviations from it close to the horizon. They exist 
for continuous small deformations 5u around the bi-Schwarzschild-AdS solution. An 
example of such solutions is shown in Fig. [l} 



ln(r/r h ) 



Figure 1. Left: Example of metric functions for an hairy asymptotically AdS solution 
for = —0.1, «4 = —0.3, and n = 0.41, where IVo, Q o, To and «o correspond to the 
Schwarzschild-AdS solution with u = 2.6333. Right: Function U'(r) for different values 
of u. Taken from [5jJ. 


Interestingly, when one takes the limit where the horizon radius vanish ry t —> 0, 
solutions with no horizon and purely made of massive field modes still exist [20] . These 
are globally regular solutions, including at r = 0, and asymptotically approach AdS. 
More recently, in addition to these solutions, Ref. [56] showed that for some discrete 
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values of u there also exist wormhole solutions which are asymptotically AdS. The 
existence of such solutions is related to the fact that the energy-momentum tensors 0 
and (12) do not, in general, satisfy the null energy condition [58]. 


3.2.2. Asymptotically flat hairy solutions. Unlike the AdS case, constructing 
asymptotically hairy flat BHs turns out to be much more complicated. In fact, as 
was pointed out in Ref. ra, solutions which differ from Schwarzschild can only exist 
for discrete values of u (fixing all the other parameters), unlike in AdS where solutions 
can be found varying u continuously. 

Those were found in Ref. [23], with the solutions having an explicit “Yukawa” 
behavior when r —>■ oo: 


N — 1 ° l + C ”-( 1 + r ^ e -v 

2 r 2 r 

(54) 

Y 1 Cl ^(l + r h) c -ru 

2 r 2 r 

(55) 

C 2 (l+r/i + rV) 

^ ,, 2^2 e ’ 

(56) 


where C\ and C 2 are integration constants. Note that the Yukawa behavior at infinity 
justifies identifying the parameter p with the graviton mass. The value of u for these 
solutions can be found fixing the values of /1, a 3 , a 4 and ft, integrating numerically from 
the horizon with the boundary conditions (52R(53), and then shooting for the value of 
u such that the solution matches the asymptotic behavior (54)—(56). Mathematica 
notebooks to generate these hairy solutions are available online EZ|. 

A trivial solution for any value of /r, 03 and ct 4 is obtained when u — 1 , and it 
corresponds to the two metrics being equal and described by the Schwarzschild solution. 
On the other hand, for 11 ^ 1 there are also regular, asymptotically flat BHs endowed 
with a nontrivial massive graviton hair. The most important result of these studies is 
that hairy solutions exist near the threshold pMs < 0.438 for any value of 0 % aq Q As 
we will discuss in Sec. 4.1 this is related to an instability of the bidiagonal Schwarschild 
solution found at the linear level for massive graviton masses satisfying precisely this 
bound [28, [29]. The existence of this instability was in fact what prompted the authors 
of Ref. [[23] to search for hairy solutions. At the threshold ~ 0.438 the hairy 
BH solution merges with the Schwarzschild solution, and above it the only bidiagonal 
solution seems to be the bi-Schwarschild one, which is consistent with the fact that this 
solution is linearly stable in this regime (see Sec 4.1 for more details). Examples of 
solutions for different choices of a 3 and ct 4 are shown in Fig. [2] 

The behavior at smaller pMs is more convoluted as it depends strongly on the 
nonlinear terms of the potential ([ 8 ]). In fact, the solutions were found to stop to exist 
below a parameter dependent cutoff p c Ms, and thus for some choices of the parameters 


+ The parameter u = U(rh)/rh remains invariant under rescaling transformations. This can used to 
express all dimensionful quantities in terms of the mass of a Schwarzschild BH with horizon r>, , i.e. 
M s = r h /2. 
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«3 and a 4 , hairy BHs with arbitrarily small fxMs seem to be excluded. For graviton 
masses /i of the order of the Hubble scale this would mean that for these parameters only 
cosmologically large BHs exist, thus unlikely to be astrophysically relevant. Through 
a detailed numerical analysis of the two-dimensional parameter space ( 03 , 0 : 4 ), these 
solutions were conjectured to follow the classification summarized in Fig. [3j However, 
it is unclear if for some parameters the solutions cease to exist when /j,Ms —> 0 or 
if they simply become extremely difficult to find numerically. In fact, in addition to 
the solution with asymptotic behavior (54)—(56), there is always another branch which 
diverges exponentially at spatial infinity. Thus, any small deviation from a regular 
solution leads to a singular behavior, making it numerically challenging to shoot for 
the correct solution when f-iMg —> 0. It might be possible that using other numerical 
methods (such as a relaxation method) turns out to be more efficient. Therefore, it is 
unclear if the classification of Fig. [3] is correct and further study is needed. 



Figure 2. Example of solutions for different values of the mass coupling y,Mg. The 
behavior is similar for any value 0:3 and 0:4 near the threshold yMg ~ 0.438 but for 
small yMg it can be very different depending on the specific values of the parameters. 
Left panel: 0:3 = —2 , a 4 = 2 and k = 1. Center panel: 0:3 = — 1, 0:4 = 0 and k = 1. 
Right panel: Metric functions for yMg = 0.01, <23 = — 2 , 0:4 = 2 and k = 1. Taken 
from [231 • 


Finally, the above picture was found to break down in the limit where one of the 
metrics is taken to be a non-dynamical Schwarzschild metric (k 1). In this case the 
numerical search suggest that, for any choice of 0:3 and 0 : 4 , hairy BH solutions exist near 
the threshold /iMg < 0.438 but they do not exist for arbitrarily small p-Mg. 

To close this Section let us point out that asymptotically de-Sitter hairy BHs 
have, so far, not been found. The study of Ref. ED] seems to indicate that any small 
deviation from the Schwarzschild-de Sitter solution develops a curvature singularity 
at finite proper distance from the event horizon. However, whether a discrete set of 
solutions exist, such as the ones discussed in this Section, is still unknown. As in 
the asymptotically flat solutions discussed above, if any small deviation from a regular 
solution leads to a singular solution, those would be very difficult to find using a shooting 
method. In this case, other methods might be more appropriate. 

3.2.3. Other special black hole solutions. Other less interesting and pathological BH 
solutions were also found in Ref. [2D]. For a matter of completeness let us give here a 
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Figure 3. Conjectured diagram of the parameter space for BHs with massive graviton 
hair in bimetric massive gravity, (i) 0:3 A — a^Ua^ = — 0.4 < — 1 - The solutions stop to 
exist below a cutoff y c Ms (which depends on the parameters); (ii) 1 < «3 = —a 4 < —2 
The solutions disappear only near pMs ~ 0.01 and are singular at small pMg, 
because some component of the metric is vanishing where the metric g^ is 
regular; (iii) 0:3 = — > —2 - The solutions exist for arbitrarily small jiMs and 

are nonsingular. Taken from 7D . 


brief description of these solutions. 


Asymptotically U, a black holes. A class of bidiagonal solutions exist where the metric 


functions U and a of the metric /, given by Eq. (47), asymptotically go to a constant 


value, while the metric g asymptotically resembles AdS at leading order in 1/r (but 
not at all orders). These solutions are generically obtained if one continuously varies 


u in the boundary condition (53) around the Schwarzschild solution u = 1, and are in 


fact asymptotically degenerate since det(/) —> 0 when r —» 00 and thus not physically 
relevant. 


Tachyonic black holes. A class of special solutions with unusual asymptotic behavior 
was also found in the massive gravity limit where the non-dynamical metric is taken to 
be Schwarzschild. In this case, solutions with the asymptotic behavior (54)—(56) exist 
which have an imaginary graviton mass /i. This suggest that around these solutions 
linearized perturbations are tachyonic, signaling an instability. 


4. Perturbations of black holes in massive (bi)-gravity 

Besides the existence of BH solutions in massive gravity, a fundamental issue one 
should raise concerns the stability properties of these solutions. Perturbation theory 
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allow us to understand these properties and can sometimes even predict the existence 
of new solutions. Recently, the first steps of this programme have been done in 
Refs. [28, 291 S3, 55], 00]. They studied linear perturbations of some of the BH solutions 
discussed in the previous sections. 

Surprisingly, they showed that the bidiagonal Schwarzschild BH solution is unstable. 
This instability was first discovered in Ref. [28] who showed that the mass term for a 
massive spin-2 held in this background is equivalent to the Kaluza-Klein momentum of a 
four-dimensional Schwarzschild BH extended into a hat higher dimensional spacetime. 
These black strings, solutions of GR, have been shown to be unstable against long- 
wavelength perturbations, the so-called Gregory-Lahannne instability [59]. This led 
the authors of Ref. [28] to conclude that perturbations of the bidiagonal Schwarzschild 
BH would generically give rise to a spherically symmetric instability. This was fully 
confirmed and extended to the bidiagonal Schwarzschild-de Sitter solution in Ref. 


Perturbations of the bi-Kerr BH solution (37) and (38), for f /JU = g^ u , were also 
studied in Ref. [29], where it was found that due to the superradiant scattering of 
bosonic helds with spinning BHs another kind of instability could be triggered. Due to 
the dissipative nature of the BH horizon and to the existence of negative-energy states 
in the ergoregion of a spinning BH, low-frequency to monochromatic bosonic waves 
scattered off rotating BHs are amplified whenever the following condition is met 


U! < mVt 


H : 


(57) 


where VIh is the angular velocity of the BH horizon and m is an integer characterizing 
the azimuthal dependence of the wave. The extra energy deposited in the wavepacket’s 
amplitude is extracted from the BH, which spins down in the process. If a confining 
mechanism is able to trap superradiant modes near the BH, this can trigger an instability 
(see e.g. Refs. [35[ i33] and references therein). It turns out that a massive bosonic held 
can naturally confine low-frequency radiation due to the Yukawa-like suppression of the 
held at large distances, ~ e~^ r jr. This instability was explicitly shown to occur for 
scalar (spin-0) [36, [23 and vector (spin-1) [38, 39], 02] helds and Ref. [29] showed that 
massive tensor (spin-2) helds would also be unstable. 

In this section we will consider small metric perturbations of the form g [LV = 
gjiJ + h where gj! J J denotes the background metric at zeroth order in perturbation 
theory and h^ denotes a small perturbation of the background metric. To study spin- 
2 perturbations h^ of a generic spherically symmetric spacetime one can decompose 
the perturbations in tensor harmonics and write them in Fourier space as (in spherical 
coordinates) [61] 

/ +oo 

e-“ T, M) + T. 9, 0)] doj . (58) 

L,m -°° 

where h a ^ a1,lm and hP° lar,Zm are explicitly given by 
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/ 0 0 h l ™{r) csc0<9^Y; m (0,0) -h l ^{r) s\n9d e Y lm {0, 0) \ 


* 0 h l ™{r ) esc 9d^Yi m {9, 0) 

* * —ho m (r) Xl m(s,<t>) 

y * * * 

C m K r »W) = 

/ H l 0 m (r)Y lm H[ m (r)Y lm 
r)Y lm 


-^^(r) sin 9d g Yi m (9,(j)) 
h l ™{r) sin 9Wi m (9, 4>) 
h l ^{r) sin 9X lrn {6,(j)) / 


(59) 


ri l o m (r)d 0 Yi 


H l ™( 


rj[ m (r)d e Yi 


lm 

lm 


V0 m ( r )d<t> Y lm 

v[ m (r)d^Yi m 


\ 


\ 


* 

* 


r 2 h t 


ee 


r 2 G lm (r)X, 


lm 


* r 2 sin 2 dh H y 

where asterisks represent symmetric components, hgg = K lm (r)Y], 


(60) 


Um + G lm (r)W lm , 

= K lm (r)Yi m — G lm (r)Wi m , Yi m = Yi m (9,(ft) are the scalar spherical harmonics 

and 

Xim(0,4>) = 2d 0 [d e Y lm - cot 9Yi m ] , (61) 

W lm {9, 0) = d 2 e Y lm - cot 9d e Y lm - esc 2 9d\Y lm . (62) 

The perturbation variables are classified as “polar” or “axial” depending on how they 
transform under parity inversion [9 —y ti — 9, (j) —> </> + n). Polar perturbations are 
multiplied by (—1)* whereas axial perturbations pick up the opposite sign (—l) i+1 (see 
e.g. Ref. [62] for further terminology used in the literature). 

In general the angular and radial parts of the perturbation h^ will fully decouple, 
such that the radial components will satisfy a set ODE’s, which together with 
suitable boundary conditions at the BH horizon and at spatial infinity define an 
eigenvalue problem for the frequency oj. These boundary conditions impose that the 
eigenfrequencies are generically complex, u j = ujr + icoj [62]. Through Eq. (58), an 


instability corresponds to an eigenfrequency with ujj > 0 with an instability time scale 
r = 1 /loj, while the case u>i < 0 corresponds to stable modes that decay exponentially 
in time. 


J h l. Radial perturbations 


Let us first consider radial perturbations of the spherically symmetric solutions (16), 


(17), and for the sake of simplicity we will mostly consider the case of asymptotically 


(63) 

(64) 


flat spacetimes, l g — If —> oo, 

ds 2 = — ^1-- j dv 2 + 2 dvdr + r 2 dQ 2 , 

ds 2 = C 2 — ^1 —— j dv 2 + 2 dvdr + r 2 dfl 2 

We will consider both the bidiagonal and non-bidiagonal solutions. The bidiagonal 
case is realized for r g = rf and either C = 1, A g = Af = 0 (the case of identical 
metrics, f tlu = g^)] or f gi/ = C 2 g fiu and with the following relation between C and the 
parameters of the action ([7]) (cf. Sec. 3.1), 


- A g , C nA f - A f . 


(65) 
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As we showed in Sec. 3.1, there is another class of solutions following from the 


ansatz (63), (64), the non-bidiagonal class, r g ^ rf. For this class of solutions the 


scale factor C satisfies the relation (26), and the parameters of the action must be 


related via the same relation as for the bidiagonal case for 1, Eq. (65). 


The metric perturbations /tjfj and , corresponding to the g and / metrics, satisfy 
the linearized field equations 

' 2 'V=9. 


5G **„ = m z 5T 


K 




( 66 ) 


For spherically symmetric modes only the terms proportional to Yi m in the 
ansatz (58) remain. Writing the perturbations in ingoing Eddington-Finkelstein 


coordinates (and changing the notation of the radial functions to avoid confusion) we 
then have 

KZJr) 0 0 \ 


Z = 


%(’•) 

h vr 


0 

0 


h r Z(r) 


( h v ,t\[r) 


IZ V - 
a (f) “ 


c 2 


V 


w f 

h vr ( 

0 
0 


‘GO' 

rr i 

\g) ( 

0 

0 


w 

h lh(r) 

0 
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0 

0 

h tt)w 

r 2 

0 

0 

0 


o 
o 

0 

r 2 sin 2 9 

0 

0 

0 

K VI 

r 2 sin 2 0 


(67) 


\ 


/ 


( 68 ) 


where an overall 1/C 2 factor for /-metric is introduced for convenience. Note that 
the advanced time v is regular at the future horizon. Therefore we require the metric 
perturbations ^(r) to be regular at r = r g and r = rf. 


4-1.1. Non-bidiagonal solutions. The calculation of the mass term in Eq. ( 66 ) yields 
the remarkably simple expression, 

/ 0 0 0 0 \ 

0 


r™ _ A ( r 9~ r f ) p -iuv 

4 r 


hf, 0 0 

hr, 

0 0 A=) 0 

hr. 

0 0 0 


s 1] = -«n, 


(69) 


/ 


where, 

A = 

and we defined, 


C 2 (j3(C - l ) 2 - 1) 


C- 1 


(70) 


"?-> - h Z - c Xr 


Notice that for the above definition, e.g. h/\(r) = hZ(r) — h/JJr), taking into account 
the factor 1/C 2 in the definition of h^y In the bidiagonal case the calculations are much 
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simpler (we consider this case in detail below), since the perturbed mass term has the 
Fierz-Pauli-like form, 


8T» V = j (<3(C - l) 2 - 2 a(C - 1) + l) e~ iuv (8^ - h^) , 


5 


= -8T» 


(71) 

(72) 


At the intersection of the two branches of solutions, for r g — rf in (69) and for C fixed 


by (26) in (71), we find ST ^ = 0. This means that the perturbation equations are those 


of GR, as it is evident from Eqs. (66). In the case of different radii, r g ^ rf, and, in 
particular, for a flat metric / the perturbations are different from GR, unless *4. = 0 
in Eq. (69). The condition A = 0, which reduces to /3 = (C — 1)~ 2 implies from (26) 
(3 = a 2 . This particular choice of parameters, which was also studied in Ref. [55], is not 


generic and corresponds to the enhanced symmetry BHs discussed in Sec. 3.1 


The divergence of Eq. (66) leads to the constraint, 


w (^f n 


°c = 0. 


Applied to the non-bidiagonal case, this constraint gives, 
A [r g - r f ) iu} 

4r 2 


K-)) ,^,0,0 =0, 


(73) 


(74) 


which immediately leads to the following simple condition on the components of the 
metric perturbations, 


h vv — n h ee — — 
n (~) — u > n (~) — > 


Co 


‘(-) 


(-) 


(75) 


where Cq is an integration constant. Using (75) back in (69) one can see that there is 


CX h^y 


only one nontrivial component of the matrix ST ,J V , namely, ST r 

In the non bidiagonal case there is no simple way to separate the gauge-invariant 
metric components from the GR part. However, thanks to the extremely simple form 


of (69) and the constraints (75) it is possible to obtain an analytical solution of the 


perturbation equations (in contrast to the bidiagonal case, where the equation (83) 
must be solved numerically, as discussed below). Note that once (75) is substituted 


into (66), we get linear differential equations (equivalent to those of GR), with the 


r.h.s. being source terms, given by the constraints (75). Therefore, general solutions for 
h'l^ and contain a part that is gauge-dependent (same as in GR) plus a particular 
solution to the full equations, i.e. 


r \g,f) ~ n GR ' (m) • 


(76) 


The particular (gauge-invariant) solution is given by a single nonzero component for 
each metric perturbation 


h rr(g) _ 
n {m) ~ 


u rr (f) _ 

n (m) 


A r g - r f) e ~ 


Aiu: 


~m 2 h e y_y 


(77) 


h 


rr(g) 

(™) 

K 


( 78 ) 
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Both the homogeneous parts of the solution, and are, individually, pure 

gauge. They can be written as Ii'q R = where to recover (67), we set 

= e~ lU}V {^°(r), ^(r), 0, 0}. The relation between and is fixed by the constraint 
(75), namely, + c i, £(/) — ^(g) + Co/2, where ci is a second integration constant. 

Since we are free to choose the gauge, we can completely eliminate the homogeneous 
(GR) part of the /-metric perturbation, i.e. 


h: 


£7 = o. 


(79) 


leaving only the non-zero , 

/ 0 


h\ 


GR 


= e 


V 


— IUJC\ 

0 

0 


— iuc\ 0 

—Co {iu + 0 

0 Cor 

0 0 


-3 


0 

0 

0 


\ 


(80) 


c 0 esc 2 (6)r 3 j 


Alternatively, one can choose the gauge such that h , Q R ' 1 = 0, but h , Q R ' > ^ 0. It is in 
general, however, not possible to set both of the homogeneous parts to zero. 

Note that static perturbations uj = 0 seem to be excluded due to the term uj ~ 1 
in (77), (78|, unless we take c 0 ~ uj. However, in the latter case, h l Q R ' 1 vanishes and 
the only nonvanishing contribution left in h rr is ~ 1/r, which describes the same non¬ 
bidiagonal solutions with properly redefined r g and ry. We can therefore exclude the 
existence of another branch of static solutions close to this family. 

Taking uj = ifl, with real positive H, one notices, that the metric perturbations are 
regular at the horizon, but they are not regular at infinity. This observation rules out 
unstable spherically symmetric modes for non-bidiagonal solutions. Moreover, even for 
real u, with “usual” boundary conditions for this type of problem, namely, outgoing 
waves at infinity and ingoing waves at the horizon, the only solution is Co = C\ = 0. On 
the other hand, non-trivial purely ingoing waves do exist (as well as purely outgoing), 
specified by two integration constants c 0 and ci, which indicates the presence of the 
hclicity-0 mode. 

Once the result in the bi-metric theory is obtained, it is not difficult to adopt it to 
the case of the original dRGT theory. For f iW = we have — o, while Eqs. (77) 
and (80) give the solution for perturbations of g /w . As in the case of the bi-metric 
theory, this solution does not have the correct behavior at r —y oo, and therefore must 
be excluded. 

Let us also comment that the scale factor C needs to be fine-tuned in the case of 


non-bidiagonal solutions by the condition (26). This particular choice of C implies the 
absence of the scalar mode for the bi-flat spacetime, i.e. the helicity-0 mode is strongly 
coupled. Indeed, from Eqs. (77), (78) and (75) one finds for r —» oo, that the “massive” 
part of the perturbations disappears. It is interesting that the scalar degree of 

freedom seems to be restored for non-zero curvature. 
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4-1.2. Bidiagonal solutions. In the bidiagonal case, the situation is very different. The 
constraint (73) requires to be traceless and divergence-free, 

= h ( _) = 0. (81) 

Taking into account Eqs. 0, 0 and the constraint fl&l| ), the perturbation equations 
(66) in the bidiagonal case can be rearranged to describe one equation for the mode 
hr *0 = + nhv^\ and another for the mode h^\ dehned earlier. The combination 

/r 0 follows the perturbed Einstein’s equation of GR, 

5G^ + k5G^ = 0, (82) 

which gives no instability, since this mode is fully equivalent to a gravitational 
perturbation of a Schwarzschild BH in GR (see e.g. [52])- On the other hand, the 
equation for /iTT is the massive Lichnerowicz equation supplemented by two constraints, 
which after generalizing it to include the cosmological constant reads: 


where A = A g = A 


□/if) + 2R otlXu h x { ^ = /i 2 /rp } , 

(/i 2 - 2A/3) h { _) = 0, 

/ and 


(83) 


y =y h + A c (p(c - 1 ) 2 - 2 a (c -1) +1) 


These equations can be shown to be consistent only if we assume the background to 
be a vacuum solution of Einstein’s equations with a cosmological constant A, so that 
R = 4A, R, IU = Ag 


jfLV 


It should be noted that under the infinitesimal coordinate 
change —y the equations of motion for the massive spin-2 perturbation 

are not invariant, while for the massless spin-2 perturbation h$ this is a gauge 
transformation, leaving the equations of motion unchanged. This is the reason why 
does not have a propagating radial mode (and hence no instability): it can be simply 
gauged away, as in GR. 

The system of Eqs. (83) has been studied by Gregory and Laflamme [59], although 


in a different context, not connected to massive gravity. They showed that this system 
admits the existence of unstable modes for 

0 ( 1 ) 


0 < /J < 


(84) 


This led the authors of Ref. [28] to the conclusion that the bi-Schwarzschild massive 
gravity solutions are unstable provided /i satisfies the condition (84). In Ref. [29] this 
instability was studied in detail also for the bi-Schwarzschild-de Sitter, and they showed 
that the unstable modes could be described by a single master equation given by 


d 2 

jpyo + [u 2 - Vo(r)] Vo = 0, 


(85) 
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where dr/dr* — f = 1 — r g /r, (po is a linear combination of polar functions in the 


decomposition (58) and the potential is given by 

1 - 2 M/r - A/3 r 2 


V n = 


x { 8 M 3 + 12M 2 r 3 (3/i 2 - 8A/3) 


( 86 ) 


r 3 [2 M + r 3 (/x 2 - 2A/3)]' 

+r 7 (p 2 - 2A/3 ) 2 [6 + r 2 (/i 2 - 2A/3)] 

— 6 Mr 4 (/i 2 - 2 A/3) [4 + r 2 (3/r 2 - 10A/3)] } , 

where we defined the Schwarzschild radius r g = rf = 2M, with M the BH’s mass (in 
units G = c = 1). This unstable mode is regular both at the horizon, where one must 
impose purely ingoing waves [62], and at infinity, where it behaves as 

(po ~ e (^>, as r —>■ oo . (87) 

The eigenvalues a; = (Ur + zo;/ for these modes can be computed by numerically 


integrating Eq. (85) with the appropriate boundary conditions. The result as a function 


of the mass coupling Mp is shown in Fig. [4} The eigenfrequencies are characterized by 
a purely imaginary (cu = icoi), positive component (loj > 0), only exist for Mp < 0.43 
and for low masses (and A = 0) behave as u>j ~ 0.7/z. 

Interestingly, for values of M and m that are phenomenologically relevant, namely 
considering the graviton mass to be of the order of the Hubble constant, m ~ H ~ 
10~ 33 eV, the mass coupling pM is always well within the instability region (assuming 
k ~ 1). This implies that a graviton with such mass would trigger an instability for 
any Schwarzschild BH with mass smaller than 10 22 M 0 ! On the other hand, for the 
same typical value of the graviton mass, as it can be seen from Fig. [4| the characteristic 
instability scale is of order of the inverse graviton mass, r ~ m -1 (and is independent 
from the BH mass). Therefore for m ~ H, k ~ 1 the characteristic instability time is 
of the order of Hubble time: for these scenarios, although the instability is present, it 
does not seem to be physically relevant for astrophysical black holes. 

The scenario of instability depicted in the previous paragraph is valid as far as the 
dimensionless parameters, in particular k, are of order of unity. However, it has been 
pointed out in [28J that the instability rate can be greatly enhanced for small k. Indeed, 


for small k, but not extremely small, so that (84) is satisfied, i.e. 


(■ r g m ) 2 < k < 1 , 

the instability rate is enhanced by a factor k~ x ^ 2 . For a given BH mass M, the maximum 
rate of instability is reached for k ~ (r 5 m) 2 , for which the instability time scale is r ~ r g . 
This means that for k ~ (M 0 m/Mj>) 2 , the instability develops in 10 -5 seconds for a 
solar mass black hole. On the other hand, in the limit k —> 0, recently considered in 
[64J for applications in cosmology, the range of instability shrinks to zero. 

As it is often the case, the existence of an instability hints at the existence of new 
equilibrium solutions branching off the instability threshold. In fact, the asymptotically 
flat hairy solutions discussed in Sec. |3.2| correspond to this new family of solutions. 
However, the end-state of this instability is still unknown, and whether it drives the 
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bidiagonal Schwarzschild solutions to these hairy solutions in some regions of the 
parameter space is unclear. 



Figure 4. Details of the instability of Schwarzschild (de Sitter) BHs against spherically 
symmetric polar modes of a massive spin-2 field. The plot shows the inverse of the 
instability timescale w/ = 1/rasa function of the graviton mass y for different values of 
the cosmological constant A = A g = Ay, including the asymptotically flat case A = 0. 
Curves are truncated when the Higuchi bound is reached y 2 = 2A/3 [65]. When 
this bound is saturated, the helicity-0 mode becomes pure gauge and the instability 
disappears [35], while below this bound the mode becomes a ghost [55] ■ For any value 
of A, unstable modes exist in the range 0 < My < 0.47, the upper bound being only 
mildly sensitive to A. Taken from [25] . 


To close, let us note that there is also a region in the parameter space for which p 2 


becomes negative. In this case, Eq. (83) does not correspond to the case studied in [59]. 
However, the negative sign of p 2 signals an instability of a more dangerous type, namely 
the spin-0 part of the graviton becomes a ghost. One can see this also from the Higuchi 
bound p 2 = 2A/3 [65J, when the de-Sitter curvature goes to zero. 


4-2. Non-radial perturbations of proportional backgrounds 

Non-radial perturbations have, so far, only been studied for proportional BH solutions. 


These perturbations are governed by the held equations (83) and were studied in detail 
in Refs. [29] 33]. Ref. [33] studied generic linear perturbations of the Schwarzschild-de 
Sitter solutions for the particular case where the Higuchi bound p 2 = 2A/3 is saturated 
and showed that the radial instability disappears, since the helicity-0 degree of freedom 
becomes pure gauge at the linear level (the so-called partially massless theory) Q 
On the other hand Ref. [[29] studied non-radial perturbations of asymptotically hat 
Schwarzschild and slow-rotating Kerr BHs, which we discuss below. 


* However, there are strong indications to believe that at the full non-linear level this symmetry is 
always broken and the helicity-0 mode reappears [66l l67l 68, 69, 70j. 
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4-2.1. Schwarzschild background. In a spherically symmetric background, when using 


the decomposition (58), perturbations with opposite parity and different harmonic index 


l decouple from each other. In addition, the radial and angular part of the field equations 
are separable such that the radial perturbation equations do not depend on the azimuthal 


number m. Inserting the ansatz (58) into the system (83) two independent systems of 


ODEs can always be found, which can be schematically written as m : 

Da*a' + V a *a' = 0, (88) 

Vp^/p l + Vp^p' = 0 , (89) 

where T>a,p are second order differential operators, Vap are matrices, and v &a, v &p 

are vectors of axial and polar functions, respectively. For l > 2, ff'A and ff'p are two 

and three-dimensional vectors, respectively. On the other hand for the dipole (/ = 1) 
mode, the angular functions Wi m and Xi m vanish and one is left with a single decoupled 
equation for the axial sector and a system of two equations for the polar sector. Finally, 
the monopole (/ = 0) mode does not exist in the axial sector since the angular part 


of the axial perturbations (59) vanishes for l = 0, while for the polar sector it can be 


reduced to Eq. (85), which, as we discussed, allows for unstable modes. The full form 


of these equations can be found in Ref. ra and is available online in a Mathematic A 

notebook I53J- 

These equations admit the generic asymptotic solution at infinity r —> oo, 




Bje 


_U r _ i 

"'OO ' j. 


+ <v 


r —> oo , 


(90) 


where koo = \Jfi 2 — a; 2 and schematically denotes the perturbation functions. This 
define two different families of physically motivated modes, which are distinguished 
according to how they behave at spatial infinity. The first family includes the standard 
quasinormal modes (QNMs), which corresponds to purely outgoing waves at infinity, i.e., 
they are defined by Bj = 0 [52] • The second family includes quasibound states, defined 
by Cj = 0. The latter correspond to modes spatially localized within the vicinity of the 
BH and that decay exponentially at spatial infinity (cf. Refs. [37, ZH [^21 SB]) - Imposing 
these boundary conditions jointly with purely ingoing waves at the horizon [62], 

$j(r) ~ e ~ lur * , r* —> -oo , (91) 


where dr /dr* — f = 1 — 2 M/r, one can numerically solve for the eigenfrequency c <j (see 
Ref. [72] for a review on numerical techniques to solve this kind of problems). In Ref. [29] 
these two different families of modes were studied in detail and no sign of instabilities 
were found, besides the spherically symmetric (l = 0) mode discussed in the previous 
section. 

Mainly due to the complexity of the polar equations, the QNMs were only studied 
for the axial sector. They found that the QNM spectrum is richer than in the GR case. 
Indeed, due to the additional degrees of freedom of the massive field, different families 
of modes exist which have no GR counterpart [|] These results are summarized in Fig. [5] 

jj For the massive bi-gravity theory, as discussed in Sec. |4.l[ there are also modes coming from the field 
equations for = hffl + nhffl , which are fully equivalent to GR. 
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M u) R 


Figure 5. QNM frequencies for axial l = 1,2 modes, for a range of field masses 
My = 0,0.04,0.52. Points with largest |wj| correspond to y —► 0. The 
fundamental mode (n = 0, circles) and the first overtones (n = 1, triangles) are 
shown. In the massless limit the “vector” modes have the same QNM frequency as the 
electromagnetic perturbations of a Schwarzschild BH in GR, and the “tensor” modes 
have the same QNM frequency as the massless gravity perturbations of a Schwarzschild 
BH. From Ref. [29|. 


Quasi-bound states were studied in both sectors and were found to follow an 
hydrogenic spectrum in the small-mass limit (M/i <C 1): 

^ / ivi a. \ 2 


u 


R 


rv.; rs_/ 




MuJr 


M/i 

/ T n T S T 1 
C sl (Mn) il+f - +2s , 


(92) 

(93) 


where n > 0 is the overtone number and S is the polarization, while the coefficient 
Csi depends on S and l. The results above are valid for moderately large couplings 
Mnv 0.2 and are in good agreement with what was previously found for massive 
spin-0 and spin-1 fields mm [38J- Note that these modes always have cui < 0 and 


thus they are stable and decay with a typical timescale r = ojJ 1 . Interestingly, Eq. (92) 
predicts a degeneracy for modes with the same value of l + n + S when M/i <C 1, which 
is akin to the degeneracy in the spectrum of the hydrogen atom. 

In addition to these modes, a new polar dipole (l = 1) mode was found [29j. This 
mode was shown to be isolated, does not follow the same small-mass behavior and does 
not have any overtone. For this mode, the real part is much smaller than the mass of 
the spin-2 held, and in the limit M/i 1 is very well fitted by 

<*W/z « 0.72(1 - M/i ), (94) 

ui/jx ~ — (M/i) 3 . (95) 


That this mode is different is not completely unexpected since in the massless limit 
it becomes unphysical. This peculiar behavior seems to be the result of a nontrivial 
coupling between the states with spin projection S — — 1 and S — 0. Besides that, this 
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mode has the largest binding energy (c Ur/ fi — 1) among all couplings M/i for massive 
fields around Schwarzschild BHs, much higher than the ground states of the scalar, 
Dirac and vector fields. 


4-2.2. Slowly rotating Kerr BHs. Massless linearized gravitational fluctuations around 
the Kerr geometry were shown to be separable in GR by Teukolsky [73] . However, the 
Teukolsky formalism does not seem to be applicable for massive spin-2 perturbations 
governed by the held equations (83). To handle the problem, Ref. [29] considered an 


expansion in the rotation parameter a = a/M = J/M 2 <C 1 where J is the BH’s 


angular momentum. Using the decomposition (58) for the metric perturbations and 


since the spacetime is not spherically symmetric, but instead axially symmetric, this 
leads to couplings between different /-modes and parity sectors, while perturbations 
with different values of m fully decouple. 

By expanding the Kerr background to first order in the spin[|J the perturbation 
equations read schematically as (cf. Ref. [38J, 72) [29]J for details) 

Ai + amAi + a{Q{Pi- 1 + Qi + {Pi + i) + 0(a 2 ) = 0 , (96) 

Vi + arnPi + a(QiAi-± + Qi+iAi+i) + 0(a 2 ) = 0 , (97) 

where, Qi = and the coefficients Ai and Vi (with various superscripts) are linear 

combinations of axial and polar perturbation variables, respectively. This system can 
always be written as a set of coupled second-order ODE’s and when a = 0, Ai and Vi 


reduces to Eqs. (88) and (89), respectively. Furthermore, it can be shown that at first 


order in the spin, couplings between different values of l do not affect the eigenspectrum 
and thus for numerical purposes one can fully separate the polar and axial sectors. 

As for the Schwarzschild case, at the horizon we must impose regular boundary 
conditions, which correspond to purely ingoing waves, 

$j(r) ~ e~ ikHr * , (98) 

as r* — >■ — oo, where = oj — mO^, Qh — a/(2Mr + ) is the horizon angular velocity 
and r + = M + a JM 2 — a 2 is the outer horizon of the Kerr geometry. All the polar 
and axial equations can be brought to a form such that the near-horizon solution is 
given by Eq. (98). If kn < 0 an observer at infinity will see waves emerging from the 
BH [73]. This corresponds exactly to the superradiance condition, Eq. (57). When 


modes satisfying this condition are confined in the vicinity of the BH, i.e. quasibound 
states, this leads to superradiant instabilities of bosonic massive fields. 

In the small M/j, limit the eigenfrequencies of the quasibound states in a 
Kerr background can be obtained from Eqs. (92)—(95) adding an extra factor 
(2 r + /i — ma/M) [21]. Thus, in the superradiant regime, the imaginary part of the 
eigenfrequency c Oj changes sign, signaling an instability. An example of these modes is 

$ As discussed in detail in [73 ESI, a second-order calculation is needed to describe the superradiant 
regime in a self-consistent way, although a first-order computation turns out to be surprisingly accurate 
in all cases explored so far. 
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presented in Fig. [6j where it is shown that the decay rate of the dipole (/ = 1) polar 
mode is very large even for small couplings M/a. Indeed, for M/i <C 1 the time scale for 



Figure 6. Absolute value of the imaginary part of the polar quasibound modes as a 
function of the BH rotation rate J/M 2 to first order in the spin for different values 
of l and m and different values of the mass coupling yM. Left panel: polar dipole 
mode for l = m = 1. Right panel: polar mode l = m = 2, S = —2. For any 
mode with m > 0, the imaginary part crosses the axis and become unstable when the 
superradiance condition ujr < mfln is met. Taken from [22j- 


this unstable mode is [29] 


T ~ 


C poUr (a/M - 2 r + u R ) ’ 


(99) 


where C vo \ ar ~ 0(1) and uj r is given by Eq. (95). This is the shortest superradiant 
instability time scale of a Kerr BH known to date. 


4-3. Superradiant instability and bounds on the graviton mass 

The superradiant instabilities discussed above have been shown to have important 
astrophysical implications, which have been recently investigated in the contexts of 
testing stringy axions and ultralight scalars [H E3 ESI CZZj, to derive bounds on light 
vector fields [39] and on massive gravitons [29] . 


Unlike the spherically symmetric instability (cf. Sec. 4.1) for which the instability 
end-state is unclear, for the superradiant instability it is clear that the BH should spin- 


down until the superradiance condition (57) is saturated m- Thus, a solid prediction 
of BH superradiant instabilities is the existence of holes in the Regge plane p3, 733, [77J 
(cf. Fig. [7J). Together with reliable spin measurements for massive BHs, this can be used 
to impose stringent constraints on the allowed mass range of ultralight bosons. These 
bounds follow from the requirement that astrophysical spinning BHs should be stable, 
in the sense that the superradiant instability time scale r should be larger than some 
observational threshold. For isolated BHs the most natural observational threshold is the 
age of the Universe, Tpubbie = 1-38 x 10 10 years. However, for supermassive BHs, possible 
spin growth due to mergers with other BHs and/or accretion also play an important role. 
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Figure 7. Contour plots in the BH Regge plane m corresponding to an instability 
time scale given by Eq. (99) shorter than a typical accretion time scale rsaipeter ~ 
4.5 x 10 7 years, for different values of the massive spin-2 field mass m g = yh (see main 
text for details). The experimental points (with error bars) refer to the supermassive 
BHs listed in [79]. Supermassive BHs lying above each of these curves would be 
unstable on an observable time scale, and therefore each point rules out a range of the 
graviton mass. Adapted from Ref. ]34] . 


The most likely mechanism to produce fastly-spinning BHs is prolonged accretion EEj. 
Therefore, a conservative assumption to estimate the astrophysical consequences of the 
instability is to compare the superradiance time scale to the minimum time scale over 
which accretion could spin up the BH. For simplicity we assume that mass growth occurs 
via accretion at the Eddington limit, i.e. assuming equilibrium between the radiation 
pressure exerted on the matter surrounding the BH and the gravitational pull of the BH, 
so that the BH mass grows exponentially with e-folding time given by the Salpeter time 
scale rsaipeter ~ 4.5 x 10' years [77]. Note that, assuming spherical symmetry, accretion 
at the Eddington rate sets an upper limit on the mass accretion rate. 

We can quantify this statement by plotting exclusion regions in the BH Regge 
plane as shown in Fig. [7] The contours correspond to an instability time scale of the 
order of the Salpeter time for four different masses of the bosonic field, considering the 
unstable mode with the largest growth rate, i.e., the polar dipole with an instability 
time scale given by Eq. (99). The plot shows that observations of supermassive BHs 
with 1O 5 M 0 < M < 1O 1O M 0 spinning above a certain threshold would exclude a wide 
range of the massive spin-2 field mass. Note that the exclusion windows extend almost 
down to J ~ 0, and this feature is important given that current spin measurements 
might be affected by large systematics. 

Nonetheless, it is clear from Fig. [7] that almost any supermassive BH spin 
measurement would exclude a considerable range of masses. Since the only parameter 
that regulates the instability is the combination fiM , similar exclusion plots exist in the 
region M 0 < M < 10 5 M o for larger values of / x Thus, the best bound comes from 
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the most massive BHs for which spin measurements are reliable, e.g. the BH candidate 
Fairall 9 [80 !|. 

Using these arguments we can obtain the following bound |29||J 

m 9 <5x 10 -23 eV, (100) 


where m g = hfi. Note that, for a single BH observation, superradiant instabilities can 
only exclude a window in the mass range of the field, as shown in Fig. [7| Nonetheless, 
by combining different BH observations in a wide range of BH masses and joining the 
superradiant bounds with bounds coming from other experiments [STj, one is able to 


constrain the range above. The constraint (100) sets a stringent bound on the mass of 
the graviton (STj, and of any massive spin-2 field governed by the field equations (83). 
If the largest known supermassive BHs with M ~ 2 x 1O 1O M 0 [E2UE3] were confirmed 
to have nonzero spin, we could get even more stringent bounds. 


5. Conclusions &; Open issues 


In this paper we reviewed BH solutions and their stability properties in massive (bi)- 
gravity. This area of research is young and fast developing — the majority of the results 
here presented have been obtained during the last few years, in the context of the recently 
discovered dRGT model. The main results include: classes of exact solutions featuring 
solutions of GR (cf. Section 3.1); numerical non-GR-like bidiagonal solutions with hair 
(Section 3.2); and perturbative (in)stability of BHs in massive (bi)-gravity (Section [4]). 
However, there are still many questions to be answered and problems to solve. We 
find it convenient to close this paper with a (necessarily incomplete and biased) list 
of outstanding open problems that are, in our opinion, still uncharted territory or not 
completely understood issues: 


As shown in ra, there is in fact an infinite family of spherically symmetric 
non-bidiagonal solutions (all related to the Schwarzschild-(A)dS metric through 
a coordinate change). However, an explicit construction of other solutions, besides 
the ones studied in Section 3.1, is still missing. More importantly, a stability 
analysis of this infinite set of solutions or a study of their global structure (along 
the lines of [50]), could help to decide which solutions are physically relevant. 

Besides the Kerr solutions discussed here, no other rotating BH solutions have been 


found so far. Extensions of BHs with massive graviton hair discussed in Sec. 3.2 


most likely exist. On the other hand the superradiant instability of the Kerr solution 
could also lead to hairy solutions analogous to the ones found in Ref. 


§ These bounds were obtained using a linearized analysis. By including the effects of gravitational- 
wave emission and gas accretion, Ref. ins showed that the linearized prediction should be corrected 
for massive scalar fields minimally coupled to GR, although such corrections would not affect the order 
of magnitude of these constraints. We expect that similar results should hold in theories of massive 
gravity. 
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• Asymptotically de-Sitter hairy solutions have not yet been found. However the 
instability of the bi-Schwarzschild-de Sitter solution is an indication that they could, 
in principle, also exist. 


• The classical no-hair theorems of GR do not trivially apply in theories with massive 
gravitons. Whether massive (bi)-gravity coupled to matter admits hairy BH 
solutions, besides the extension of standard (electro)-vacuum solutions presented 
here, is unknown. For example, using the linear Fierz-Pauli theory, Ref. [85 ] 
suggested that when the matter’s energy-momentum tensor has a non-vanishing 
trace BHs should have hair. However, it is not clear if this is still true for the full 
nonlinear theory. 

• The spherically symmetric instability uncovered in Refs. [281 29] would presumably 
cause the Schwarzschild spacetime to evolve towards another spherically symmetric 
solution. However, depending on the parameters of the theory, different types of 
solutions exist making the end-state of this instability unclear. This in turn makes 
nonlinear time evolutions of bidiagonal Schwarzschild BHs highly desirable. 


• On another related note, the non-uniqueness of the solutions in massive (bi)- 
gravity makes it unclear what is the outcome of gravitational collapse. In highly 
symmetric scenarios the bidiagonal Schwarzschild solutions seem to be the more 
natural outcome. But it is of course possible that, in some regions of parameter 
space, Schwarzschild BHs are not the preferred outcome of gravitational collapse or 
even that BHs do not form in massive gravity. These issues can only be addressed 
by performing nonlinear collapse simulations. 


• It was argued in Ref. [86] that another type of instability of black holes in massive 
gravity (with one fixed metric) should exist. Such a process would be similar to the 
discharge of a collapsing charged star, but in this case the black hole’s mass would 
decrease in a finite time (even classically). A complete numerical analysis should 
be put forward to address this question. 


The stability of the hairy solutions discussed in Sec. 3.2 remains an open issue 


• The absence (or presence) of ghosts in the perturbation’s spectrum of BH solutions 
in massive (bi)-gravity is an open issue. Although the construction of the massive 
interaction terms in dRGT theory and its extension guarantees the absence of the 
Boulware-Deser ghost, the other remaining propagating modes can possibly become 
ghostly in some backgrounds and must therefore be examined in detail. 


The linear stability of non-bidiagonal solutions to non-radial perturbations has not 
been studied. For example, whether the non-bidiagonal Kerr solution discussed in 
Sec. 3.1 suffers from a superradiant instability is unknown. 


It is possible that some of the theoretical issues discussed here might help to decide 
whether massive (bi)-gravity is a viable theory to describe our physical world or not. 
Thus, we hope that this paper will serve as a guide for future developments on this 
relatively new and exciting topic. 
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